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Abstract 

In this paper, we introduce a fictitious dynamics for describing the only fast relaxation of a stiff 
ordinary differential equation (ODE) system towards a stable low-dimensional invariant mani- 
fold in the phase-space (slow invariant manifold - SIM). As a result, the demanding problem 
of constructing SIM of any dimensions is recast into the remarkably simpler task of solving a 
properly devised ODE system by stiff numerical schemes available in the literature. In the same 
spirit, a set of equations is elaborated for local construction of the fast subspace, and possible 
initialization procedures for the above equations are discussed. The implementation to a detailed 
mechanism for combustion of hydrogen and air has been carried out, while a model with the 
exact Chapman-Enskog solution of the invariance equation is utilized as a benchmark. 
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1. Introduction 

Models for describing detailed reaction mechanisms of hydrocarbon fuels and biochem- 
ical processes in living cells are typical examples of large multiscale dynamical systems 



|25lll5|,|28|,|43|]. In this respect, modern research has to cope with an increasing difficulty mainly 
in two aspects: First, the number of degrees of freedom is tremendously large making it difficult 
to obtain a physical understanding of the above phenomena. In addition, computations are often 
dramatically time consuming due to a wide range of time-scales to be resolved. As a result, 
methodologies able to tackle the above issues become highly desirable. The issue of physical 



understanding is drawing an increasing attention in the realm of kinetic modeling of biological 



systems with many degrees of freedom |32l 1351 14611 . Modern simplification techniques (often 



referred to as model reduction methods 12011 ) are based on a systematic decoupling of the fast 
processes from the rest of the dynamics, and are typically implemented by seeking a low dimen- 
sional manifold in the phase-space. Towards this end, several methods have been suggested in the 



literature II 1 8L 1 1 QL 13 8L 1291 147L 1341 1261 1 111 which are based on the following picture (for dissipative 
systems with unique steady state to be addressed below): Multiscale systems are characterized 
by a short transient towards a stable low-dimensional manifold in the phase space, known as 
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the slow invariant manifold (SIM). The subsequent dynamics is slower and proceeds along the 
manifold, until a steady state is reached. 

Recently, the Relaxation Redistribution Method (RRM) has been proposed J5| and imple- 
mented in realistic combustion mechanisms for hydrogen (12] and methane mixtures RRM 
has been introduced as a particularly efficient scheme to implement the film equation of dynam- 
ics (see section [2] below), which can be used to construct the SIM and adaptively choose the 
minimal description of complex multiscale systems Il2l 1 1311 . In the latter References, the min- 
imal description is understood as the minimal dimension of a convergent SIM (by RRM). For 
completeness, we stress that an alternative approach for the adaptive simplification of multiscale 
systems is the G-Scheme in 14811 . 

In the present work, following the basic idea behind the RRM, we derive a set of ordinary 
differential equations which approximate the RRM dynamics (here, referred to as governing 
equations of the linearized RRM). A remarkably easy implementation of the latter method for 
constructing SIM in any dimensions is then proposed. 

This manuscript is organized as follows. In section [2] the problem of model reduction is 
introduced and the notions of both invariance equation and film equation are briefly reviewed. 
The governing equations of the linearized RRM are presented in section [3~T1 where the link 
between the presented method and other approaches (i.e. direct solution of invariance and film 
equations 1 1_8], ILDM f38ll . CSP (2^1 ) is shortly discussed. A novel algorithm for approximating 
the fast subspace is introduced in section [3T21 The accuracy in describing the SIM by governing 
equations of the linearized RRM at steady state is discussed in section [331 whereas a possible 
initialization of them is proposed in section 13.41 The suggested linearized RRM is tested in 
section|U while conclusions are drawn in section|5] 



2. Background 

Let an autonomous system of ordinary differential equations 



dY 
dt 



fi(Y) 



f„(Y) 



f(Y), 



(1) 



describe the time evolution of a state Y = [c\,...,c„] T in the phase-space U, where n is the 
dimension of U and the superscript T denotes transposition. Model reduction techniques enable 
the construction of a simplified ODE system 



dt 



= /'(£), 



(2) 



where the state £ = , £ 9 ] belongs to a reduced space S with dimension q « n, and evolves 
in time according to the slow dynamics of system ([TJ. 



2.1. Slow invariant manifolds 

By analogy with classical thermodynamics, a reduced model © represents a macroscopic 
description of a physical phenomenon (given by ([TJ) where various processes with disparate 
timescales occur. 
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Formally, the link between the microscopic world and the corresponding macroscopic de- 
scription can be established by resorting to the notion of slow invariant manifold (SIM). In other 
words, the reduced (macroscopic) dynamics (|2} occurs along a ^-dimensional SIM, Q. S,M , em- 
bedded in the phase-space U. Thus, through Q SIM , it is possible to pick up the most likely 
microstate among all the possible ones which are consistent with a macrostate characterized by 



the macroscopic observables £ e S (see also Bill ). In the following, by model reduction, we 



mainly refer to constructive methods of both the slow and fast subspaces, and we assume that an 
arbitrary manifold Q. can be defined (at least locally) by means of a mapping 

T-.E^U. (3) 

By definition, Q. S,M c U is invariant with respect to the system ([TJ if inclusion T(fo) e f2 implies 
that Y(t) € Q. for all future times t > to. Equivalently, if the tangent space T y to £2 is defined at Y, 
invariance requires: f(Y) e T y . 

In order to transform the latter condition into an equation, it proves convenient to introduce 
projector operators. For any subspace T y , let a projector P onto T y be defined with image imP = 
T y and P 2 = P. Then the invariance condition can be expressed as: 

A = (/-P)/ = 0, (4) 



where A is often called defect of invariance III 811 . and / represents the identity matrix 



It is worth stressing that, although the notion of invariance discussed above is relatively 
straightforward, slowness instead is much more delicate. We just notice that invariant mani- 
folds are not necessarily suitable for model reduction (e.g., all semi-trajectories are, by defini- 
tion, one-dimensional (ID) invariant manifolds). For singularly perturbed systems, the notion 
of slow invariant manifold has been defined in the framework of the geometric singular pertur- 
bation theory by Fenichel 11611 . However, we should also point out that in general the different 
methodologies proposed in the literature for model reduction purposes are based on different ob- 
jects. For instance, it is known that the rate controlled constrained equilibrium (RCCE) manifold 



1 271 12411 typi cally does not even fulfill the invariance condition (0), whereas other methods (see, 



pi cally 

e.g., 11181 1511. 14011 . Il38ll ) attempt the construction of invariant objects (with different accuracy). 



Here, we follow the rationale behind the Method of Invariant Manifold (MIM) HHH], where 



slowness is understood as stability (see also chapter 4 of 11811 ). so that a SIM is the stable sta- 
tionary solution of a relaxation process (film equation) 

—^■ = d-P)f. (5) 

We notice that the projector operator P introduces first order spatial derivatives (with respect to 
the manifold parameters £). Therefore, © and (0 are partial differential equations (see also 
lfl8 . 391) whose unknown is the function T, which is conveniently utilized for a parametric 



representation of the manifold CI, with Q, being the image of f: D. = T(£), ^eS. 

Several numerical schemes have been suggested in the literature for solving Eqs. J4) and 
((5): The Newton method with incomplete linearization and the relaxation method in 11811 . the 
semi-implicit scheme in olR represent a few examples. The latter approaches aiming at the 
direct solution of both the invariance condition and film equation are often hindered by severe 
numerical (Courant type) instabilities lfl4l fl9i I7I1 . Toward the end of overcoming the latter issues, 



the Relaxation Redistribution Method (RRM) has been recently introduced IU3L 11211 (see also 
Fig. [TJ. In the following, exploiting the rationale behind the RRM, we devise a set of ordinary 
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Figure 1 : (Color online). Schematic representation of the basic idea behind the Relaxation Redistribution Method (RRM). 
In a small neighborhood of the pivot Y k (large circles), a linear approximation of the SIM at the iteration k is considered. 
With the aim of driving the pivot towards the SIM, RRM I13tll3, tql prescribes an updating rale Y k — » Y k+l as schemat- 
ically sketched above (small circles represent neighbors of the pivot). Here, an ODEs system 1211 whose dynamics 
approximates the latter updating rale is suggested and tested. 



differential equations approximating the dynamics of the film equation (0 in a neighborhood of 
a fixed macrostate f eS. 



3. Linearized Relaxation Redistribution Method (RRM) 

3.1. Slow subspace 

Let the dynamical system (Q} be characterized by a hierarchy of time scales, and let t be of the 
order of the fastest scale of Let the q x n matrix B and its j'-th row Bj 



B 



b„ 



Bj = [b n ,...,b jn \, 



(6) 



define a linear mapping from the phase-space U into a reduced space H of dimension q « n: 



BY = £ 



6 



(7) 



such that a macrostate £ can be associated with any microstate Y via the (0. In the following, we 
develop an iterative methodology for refining an initial approximation of the SIM in a vicinity 
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of a given macrostate £ € S. To this end, at each iteration k, we assume that, in a neighborhood 
of f , the SIM is approximated by an affine linear mapping, T k ' : H — » [/, of the form: 



where A k and Z A are a « x g matrix and a n x 1 column vector, respectively, such that: 





4 - 












l k 

L *>j J 



(8) 



(9) 



Notice that, the over-bar denotes the pivot (Y k in Fig. [TJ at an arbitrary iteration k along with the 
corresponding macrostate 

For the state Y k = T k Q) belongs to the space defined by the linear function in ©, the column 
vector l k satisfies: 

l k = Y k -A% l = BY k . 

Assuming the existence of a ^-dimensional SIM, £l s,M , we aim at devising a procedure for 
updating the linear mapping (8): 



k+l . 



Y = A k+l % + I 



(10) 



such that dTOb describes Q s/M with a better accuracy than ((8) in a neighborhood of ^. Toward this 
end, we follow the rationale behind the Relaxation Redistribution Method (RRM) introduced in 
fl. [l2ll . We stress that, at any iteration k, the manifold is described by the mapping T k with the 
(q + 1) microstates Y , Y k + A k d£ belonging to the affine subspace ([8]). 

In Fig. [TJ we pictorially show the relaxation of Y k (large circle) and one of its q neighbors 
(small circle) 

Y k +A k d%, i=l,...,q, 

(all in the space defined by T k ) toward £l s,M during the time r, where d£ is a small parameter. 
Notice that, owing to arbitrariness in picking the ;-th neighbor, for simplicity we make the choice 
d^i = dg, Vi. According to the RRM algorithm 111211 . describes the subspace defined by the 
set of q + 1 relaxed states. The updated points can be written as: 



Y k +f(Y k )r, 

Y k +A k d£ + f(j k +A k d£)r, i=l,...,q, 



(11) 



which represent the advance in time of the (q+ 1 ) points Y , Y k +A k dg during a period r, according 
to an explicit Euler scheme. Upon linearization of the right-hand side of ([TJ, f (Y + dY) ss 
f(Y) + J (Y) dY, ( ITTb take the approximate form: 



(12) 



F* + /(f*)t, 

f* + A k d£ + f (?*) t + J (? k ) A k d%T, i = 1, q 

with J = ^(y A ) = [d/i/dy,] denoting the Jacobian matrix evaluated at Y k . Thus, a set of q 
vectors spanning the linear space described by f k+1 reads as follows: 



A k +j(f k )A k T, i=l,...,q. 
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(13) 



The RRM algorithm |12] introduces a fictitious dynamics such that an arbitrary state Y k = f k 
"moves" to a new location: 



which obeys the relationship: 



l=B j Y k = B j Y k+ \ 7 = 1, 



(14) 



(15) 



Equations ([T5T l stipulate that the movement ? A — > is orthogonal to the parameter space E 
(see also Fig. [TJ and, upon substitution of (TBI in the (Q3), enable the explicit computation of the 
variables a, in (fl4l by solving the linear system: 



1 + /J 1 /A!t 



B ci JA\t 



BiJA k r 



1 + VA*t 





«1 




Bif 






— —T 






. a q _ 




. B q f . 



(16) 



where due to ^ and (0 Z?Ay = <5y Vfc, with 5y being the Kronecker delta.In the following, 
in order to save notation, we assume that both the vector field / and the Jacobian matrix J are 
computed at Y k \ equation ( fT4b takes the more compact form: 



yk+\ =f k + (f - M k <S>l l Bf) t, 



(17) 



with denoting the q x q matrix on the left-hand side of dT6b . whereas M k is a n x q column 
matrix defined as follows: 

M k = \a\ + JA\t, ...,A k q + JA k qT ] =(I + Jr)A k , 

with I denoting the n x n unit matrix. It is worth stressing that the solution of the equations (TTtT > 
requires an updating rule for the matrix A in <j9j as well: A k — > A k+l . For this purpose, we notice 
that an arbitrary point Y in the linear space described by the function y ri+I takes the form 



r k+l : Y = ^ ai (A* + /AJt) + f * +1 , 



(18) 



hence the y-th parameter corresponding to Y is 



or, equivalently, in matrix notation: 



7 = M*A + y* +1 , 



(19) 



with A = fori, a 9 J . Equation ( fT8l stems from (TT~4-b where the origin, F* + f(Y k )r, of the affine 
subspace has been replaced with Y k+l . 
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Equations $19[ yield the function <f~ k+l in the form (TTOb : 



k+i 



where 



= M k a>- k \ i k+l = y^ 1 - M k <s>- k l l 



(20) 



The updating rules ([P71 i and ( f20T > can be interpreted as the explicit Euler numerical scheme for 
solving the following dynamical system: 



lL = f-M® l Bf, — = \m<S>- 1 -A]t- 
dt J J fit L 1 



(21) 



dt rff 

The second equation in (f2Tb can be derived, by analogy with (TTTT i. after recasting the first equation 
in (Got as follows: 

A *+i = A k + \^M k <&l x - A k ) T - l \ T. 

The above equations (fJTJ are the governing equations of the linearized RRM, which dictate a 
fictitious temporal evolution of a state Y and anxq matrix A 



[Au...,A q ], 



(defining an affine linear mapping of the form T : Y — At; + I) towards the corresponding q- 
dimensional slow invariant manifold Q S,M in a neighborhood of a given macroscopic state 
with 

M = \Ai + JAit, ...,A q + 7A 9 tJ , 0(/, j) = Sij + BjJAjT. 

For the sake of clarity, we point out that the dynamics (EH (as well as the RRM dynamics in 
1 12[]) is referred to as fictitious because, unlike the original detailed system (Q]i, no physical or 
chemical processes are typically described by it. 

Moreover, we stress that the presence of the time-scale t in the right-hand side of the equations 
in (1211 introduces a remarkable stiffness, thus the ODEs (f2TJ typically require state of the art 
stiff integrators (see, e.g., ll42"lo . According to the RRM method, the SIM is obtained when 
the relaxation and redistribution steps balance each other (details can be found in II 1210 . In the 
suggested algorithm, the analogous condition is satisfied at the steady state (here denoted as 
Y = Y ss and A = A ss ) of the dynamical system (|2TT >. Hence, the SIM is given (in a vicinity of £) 
by: 

q sim . Y = A ss £ + l ss , I s " = Y ss - A ss l (22) 

We stress that computation of the quantity /" does not require additional refinements, and it is 
performed by (l22t (upon convergence of (|2H ) if a linear approximation is to be provided for 
approximating the mapping ^ in a neighborhood of Y ss . 

It is worth noticing that, inspection of the right-hand side of the first equation in (fJTJ reveals a 
clear connection between the RRM method introduced in IU3U12I1 and the film equation (0. In 
fact, although ( |2TI ) represents a system of ordinary differential equations whereas (0 is a partial 
differential equation, the former only describes the (0 locally in a vicinity of a macrostate ^. 
In this respect, the projector onto the tangent space of a manifold Q. takes the explicit form: 
P = M<S> X B. In this respect, the latter operator satisfies the condition of projectors: 



P 2 = MOr'BMcD 



A B = M<S>~ l B = P, 



due to the relation BM = <S>. Similarly to (0, the governing equations of the linearized RRM 
(|2TT > prescribe a composition of two motions: the first one along the detailed dynamics /, while 
the second one along the tangent space of Q., -Pf = -M<f>~ l Bf. Finally, at steady state of ( |2TI ). 
the invariance condition © is satisfied: 



/ - M<S>~ l Bf = (l- M<5 _1 b)/ = (/ - P)f = 0. 



The above equation imposes that, on the SIM, the component of dynamics / in the fast subspace 
vanishes. Since that condition lies at the heart of other popular methods (such as ILDM and CSP 
1 38, 2910 . this explains the formal resemblance of (t2TT > (at steady state) to the equations adopted 



in ILDM and CSP. 



3.2. Fast subspace 

The methodology proposed in the previous section can be utilized for extracting the slow 
invariant manifold (i.e. the subspace of slow motions or slow subspace, for short) with respect 
to the ODE system ([TJ. Nevertheless, this is only one aspect of model reduction: Computing 
the fast subspace is indeed required in order to achieve the complete decomposition of the full 
dynamics / (slow-fast decomposition). 

We notice that, towards this end, several approaches have been proposed in the literature. For 
instance, the notion of thermodynamic projector for dissipative systems supporte d by 

a thermodynamic Lyapunov function, the spectral decomposition of the Jacobian matrix J 13811 . 
and the CSP algorithm l22ll are some of the most popular examples. 

Those methods might be adopted in combination with the above technique (l2TT i as well, for an 
a posteriori reconstruction of the fast subspace. However, here in the same spirit of the method 
presented in section |3~T1 we propose an alternative procedure for computing the fast subspace, in 
a neighborhood of a given macrostate once the linear function (l22l has been computed. 

Let us assume that the fast subspace can be uniquely parameterized (at least locally) by the 
variables 77,-, i = 1, n — q — r. Let the z x n matrix B and its j-th row Bj 



[bjU-,bjn]. 



(23) 



define a linear mapping: 



BY — 77, 77 = 



(24) 



where the dynamics / of the system ([T) obeys a set of r linear conservation laws. In a neighbor- 
hood of the SIM point Y ss , at a given iteration k, the fast sub-space can be represented by a linear 
function as follows: 

Y = A k ij + f, (25) 



with 





r ~a k ■ 
"11 






2* 




n- 


A k = 










, f = 






a k . ■ 

L ni 


■ 4- 




a k . 

L ni J 







(26) 




Figure 2: (Color online). Rationale behind the refinement process of the fast subspace: In a neighborhood of the SIM 
(slow subspace), the anti-parallel dynamics — / reacts with a torque if (25) does not span the fast subspace. As a result, 
the latter subspace is the stable stationary solution of Eqs. i30\ and )3U . Small circles denote neighbors of a pivot (large 
circle). 
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Similarly to the procedure of section |3~T1 here we aim at devising an iterative procedure 

so that the function d25l l. in the limit k — > oo, accurately describes the fast subspace. Let drj 
denote an arbitrary small parameter. Following the pictorial representation of Fig. [2] for every 
variable 77, we consider the relaxation of the two neighbors (in the affine space d25l )) of Y ss 

Y ss + A*di], Y ss - A k drj, 

under the anti-parallel dynamics -/. After time r, these states move to the new locations: 

Y ss + A k drj - / (Y ss + A k djj) t, Y ss - A k djj - / (Y ss - A*drj) t, 

which, upon linearization of the vector field /, take the approximate form: 

Y ss +A k i dT]-f(Y ss )T-J(Y ss )A k i dr r r, Y ss - A k dt] - f (Y ss )t + J {Y ss )A k i dr r r. 

Therefore, a set of z vectors spanning the fast sub-space at the iteration k + 1 reads: 

A k -J ss A k T, i=l,...,z, 

where, for the sake of notations, J ss = J(Y SS ). We can thus describe the updated fast sub-space 
as follows: 

Y - Y ss = J] &t (if - J ss A k r) , 

i=l 

or equivalently in matrix notations 

Y - Y ss = M k A, (27) 

with 

M k = [A\ - J ss A\t, ...,A k - 7 ss A*tJ , A = [a u ...,af . (28) 
By substituting the (|27) into d24b , 

Tj - BY SS = <J> A A, (29) 
where the generic element 3>t (;, J) of the zx z matrix reads 

O k {i,J)=Stj-BiJ"A)T, 

where, owing to d24l i and d25l l. B,A* = dij Vk. From d27] > and d29l , it follows that at the iteration 
k + 1 the linear function describing the fast sub-space is: 

Y = M&l 1 (77 - BY SS ) + Y ss , 

so that 

A k+l =M k % x . (30) 

Similarly to the (f20b . the updating rule (f3Qb can be regarded as the explicit Euler scheme for 
integrating the dynamical system: 

— =\m®- 1 -A]t-\ (31) 
dt L J 

with 

M = [Ii - / ss Air, A, - J ss A z t] , O (i, j) = dij - BjJ ss AjT. 
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3.3. Stationary solutions of the linearized RRM equations 

We notice that the matrix A keeps evolving under the dynamics (l2TT i until the following steady 
condition holds: 

dA 

— =M0 _1 -A = O, (32) 
dt 

which can be recast in the more explicit form: 

(/ + Jt)A(I + BJAtY 1 = A. (33) 

It is straightforward to prove that right eigenvectors of the Jacobian J satisfy the stationary con- 
dition ( l33l . Let the columns of A represent a set of q eigenvectors of / such that 

J A = AL, (34) 

where L is a qxq diagonal matrix whose non-zero components are the corresponding eigenvalues. 
Upon substitution of ( |34| ) in ( l33l . we obtain the identity: 

A(I + Lt)(I + BALtY 1 = A, (35) 

due to the condition BA = I. Similar considerations apply to the evolution of A under the 
d3TT >. Hence, we can conclude that the eigenvectors of J (evaluated at Y ss ) do provide stationary 
solution for both the equations dA/dt = and dA/dt = 0. The identity d35l ) also suggests that, 
if ( 1341 ) holds, the projector operator in the first equation of OTT) takes the simple stationary form 
P = MO'B = AB, such that the pivot evolution is ruled by: 

% = (I-P)f = {I -AB)f. (36) 
dt 

The above considerations suggest that the proposed method can deliver approximations of the 
SIM up to an accuracy of the order of ILDM [38]. 

It is worth stressing that such a limit is not due to the RRM approach 

II EI El, 

rather to 

the linear approximations of (|3) and ([TO by (|8) and (fTZt . respectively. Hence, other governing 
equations leading to more accurate description of the SIM compared to (f2Tb may be also devised, 
abandoning the present linear expressions (|8} and (fTZt in favor of higher order approximations 
(at the cost of a more demanding implementation). 

3.4. Initialization 

The method described in section IXTl for constructing local approximations of SIM requires the 
initial choice of Y l and A 1 . Several approximations of the SIM can be adopted for this purpose as 
proposed in 1 3(1 33 ] . In the following, we discuss in detail another possible initialization strategy 



for the case of dissipative systems. Closed chemically reactive mixture of gases are prototypi- 
cal examples of large dissipative systems that can be addressed by model reduction techniques 



1 2111 . In fact, due to the second law of thermodynamics, the dynamical system (Q3, describing 
the temporal evolution of chemical species concentrations, is equipped with a thermodynamic 
Lyapunov function related to entropy and always decreasing in time. In this case, a rough 



approximation of the SIM is often provided by the quasi-equilibrium manifold (QEM) 1 18, 71 [ 
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also referred to as constrained equilibrium manifold I27I B Pioll . A QEM is defined by means of 
the following constrained optimization problem: 



' G (F) -> min, 

• BF=[|i,...,? ? ] r , (37) 
. DY=\xu..,Xr] T . 

where q denotes the QEM dimension, while the r x n matrix D imposes the conservation of 
the number of moles (Xi,i=i,...,r) of r chemical elements involved in the reaction. In Fig. [3] 
the geometry behind the notion of QEM is shown schematically. Let H and B be the second 
derivative matrix of the Lyapunov function G and the null space of the (q + r) x n matrix of 
contraints in (l37l i: 

Di 



Let the n x (n - r) matrix D be defined as follows: 

D = [D u . ..,£>„ - r ] = kerD. 

With the basis vectors \D\, D„-ri spanning the null space of D (kerD), an arbitrary vector 8Y 
along the tangent space of a QEM can be written in terms of the vector 5 = [5\ , 6„- r ] T as: 



B = ker 



#1 



Di 



D r 



D = 



D6 = 6Y = D 



Sn-r 



(38) 



The geometry behind the optimization problem (l37t imposes the orthogonality condition 
B T HD5 = (see also Fig. [3]), and the tangent space to a QEM is spanned by: 

T QEM =Dker(B T HD). 

Recalling the definition of the matrix A in (O, a possible initialization of A (with T l : Y — 
A l £ + /' describing locally a QEM) takes the form: 

A 1 = T qem {BT Q emT\ (39) 



whereas f 1 J = Y l can be found by solving the optimization problem (l37l i using for example 
tools suggested in 1 41 ] and J^]. Finally, possible choices for the matrix B are discussed in lUQl 
(spectral quasi equilibrium parameterization) and 9271 12[ (constrained equilibrium parameteriza- 
tion) while exact formulae for computing matrices H and J can be found in J5, LI]. 

Moreover, Eq. (l30t and the dynamical system (l3Lt require the initial condition A 1 = A(t„ = 0). 
A possible option is the following: Since fast motions are necessarily transversal to the slow 

12 



G iso- 




—k ► & 

Figure 3: (Color online). Pictorial representation of the notion of quasi-equilibrium manifold (QEM) )37t . H and B 
being the second derivative matrix of the Lyapunov function G and the null space of the full set of constraints in )37t . 
respectively. 
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subspace d22l i. a reasonable choice for the (n — q — r)X n matrix B reads: 



B 



ker 



(40) 



D, 

As a first guess (k = 1), let the mapping d25b describe the orthogonal subspace to the SIM d22b . 
More specifically, let the former space be spanned by the columns of the nX(n - q - r) matrix 



T- 1 =Z)ker[(A ss ) 7 D], 
similarly to ( ]39l l, initial condition for (]30l and (|3TT i takes the explicit form 

A 1 = T ± [BT ± )~ 1 . 



(41) 



(42) 



We notice that owing to the relations BA 1 = / and BA l = I, as an alternative to (f39l > and d42l . 
A and A can be initialized by computing the Moore-Penrose pseudoinverse matrices of B and B, 
respectively. 

3.5. Stability and adaptive construction of slow invariant manifolds 

Stability of the governing equations of the linearized RRM d2"TT i can be exploited for adaptive 
construction of SIM. In the first place, Eqs. (l2"TT i can be solved with q — 1: If convergence 
is experienced, we assume that a ID reduced model of the system ([TJ can be constructed in a 
vicinity of the macrostate In other words, a minimal description (f2]i of the detailed system ([TJ 
can be accomplished by means of one degree of freedom. On the contrary, with no convergence, 
the manifold dimension is updated to q = 2 and the procedure repeated. Upon convergence 
with some q = q, we may infer that a minimal description of the detailed dynamics requires q 
degrees of freedom. In this sense, the suggested method enables an adaptive construction of SIM 
(i.e. varying dimension in the phase-space without any a priori assumptions on the value of q). 
The above idea relies upon the assumption that RRM is stable provided the existence of SIM of 
a certain dimension q 113, 12]. More details on the stability of the RRM can be found in J5|, 
where a comparative study between a method for the direct solution of the film equation (01 and 
RRM is performed. 



4. Benchmark 



For the sake of simplicity, we consider here a four-dimensional model where the dynamics of 
two fast variables (C3 and C4) is slaved to the motion of the slow variables c\ and C2 B4411 . Let 
the functions f, f%, Q\ and 82 depend on cj and ci only. In the following, we focus on the ODE 
system: 

f (ci,c 2 ) 

dc fi{c\,c 2 ) 
dt ~ -\ [c 3 - 0i (ci,c 2 )] +fid Cl 6\ (ci,c 2 ) +fzd C2 0i (ci,c 2 ) 
-- [C4 - 02 (ci, c 2 )] + fd Cl 6 2 (ci,c 2 ) + f 2 d C2 6 2 {c\, c 2 ) 
14 



(43) 



□ INITIAL PIVOT 




Figure 4: (Color online). Slow invariant manifold with respect to the dynamical system (43) with (49), u> = 3 and 
e = 0.025. Starting from the initial conditions i50i , the governing equations of the linearized RRM 1211 with n = 4, 
q = 2 and r = 3 X 10~'° are solved by means of the stiff numerical scheme odel5s readily available in Matlab® l42ll . 
Here, the steady state is reached after an integration time Tf = 1, At steady state, the solution trajectory (dots) finally 
lands on the SIM. 




Figure 5: (Color online). An array of initial states has been refined by means of the linearized RRM. Stationary states of 
the d — 1 1 are reported (circles). For a comparison, the exact Chapman-Enskog solution and a detailed solution trajectory 
of the system (43) are also shown. Here we use ca = 3, r = 3 X 10~ 10 , while integration of 12 U is performed for an 
integration time Tr = 1 at any point. The computational time required to refine the entire array, composed by (21 X 21) 
points, was 2.5 minutes using a Matlab code on a single processor with 1.73 GHz. 
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Deviation of c 




Figure 6: (Color online). Difference between the SIM constructed by the linearized RRM method as reported in Fig. [5] 
and the Chapman-Enskog solution (C3 = 6\, c 4 = 9j). 



where e and 3, denote a fixed small quantity and partial derivative with respect to variable i, 
respectively. Assuming that the dynamics of c 3 and c 4 is slaved to the slow variables, 



c 3 = c 3 {c\, e 2 ) , c 4 = c 4 (a , c 2 ) ■ 
and, according to the chain rule, time derivatives of (l44t take the explicit form: 

% = ^ Cl C3^ + 3 C2 C 3 ^ = d fl C 3 /l + 5 C2 c 3 / 2) 



~3T = d ci c *W + d a c ^iiT = ^ci c 4/l + d C2 C 4 f2- 



(44) 



(45) 



Upon substituting equations d45l > in ( l43l l. one obtains the following invariance conditions with 
respect to (l43l l (see also section|2]and Eq. (0]i): 



i [c 3 -0i (ci,c 2 )] +fid Cl 8i (ci,c 2 )+/ 2 3 C2 0i (ci,c 2 ), 



<9 Cl C 3 /l +<9 f2 C 3 / 2 

<9 Cl c 4 /i + 3 c ,c 4 / 2 = -i [c 4 - 2 (ci, c 2 )] + /i3 Cl 2 (ci, c 2 ) + fid Ci e 2 {c\, c 2 ) 



A common approach to obtain solutions to the above invariance conditions is the Chapman- 
Enskog method 0,1451], which is based on the assumption that e is small compared to all other 
quantities, and it is implemented by series expansions of the (l44l i in powers of e: 



c 3 ( Cl ,c 2 ) = cf + ec* 1 ' + e 2 c< 2) + .. 
c 4 (c u c 2 ) = cf + ec < 1) + e 2 4 2) + .. 



(47) 



Hence, the first equation in (|46| | reads: 

M, [cf + + e 2 cf + ...] + f 2 d C2 [cf + + e 2 ^ + ...] = 
_I [4°> + ec™ + e 2 4 2) + ... - 9i] + /i3 Cl 0j + / 2 <9 C2 0i. 
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Figure 7: (Color online). On the basis of the system J43t with I49K w = 3 and e = 0.025, we solve )2U by means of the 
numerical scheme odel5s |42T| starting from several initial conditions Y l and A 1 , with fixed parameters f i = ci = 0.3 
and f2 = C2 = —0.2. Time evolution of C3 and C4, as dictated by 42H . is reported. For an integration time Tt > 1, the 



steady state is reached such that c\ 



-0.4422, c" = 0.2520 for any initial conditions. At steady state, small deviations 



from the Chapman-Enskog solution (star) are observed (smaller than 0.01%). 



After collecting terms with the same power of e, we obtain: 



h (ci,c 2 ), 



0, Vi > 0, 



(48) 



namely C3 = f?i (c\,C2), and similarly C4 = 62 (ci, C2). 
For illustration purposes, we choose: 

01 (ci,c 2 ) = sin(wci)sin(wc 2 ), 02(ci,c 2 ) = [(1 +e"" C| )(l + e" 
/1 = -ci, fi = -2c 2 . 



(49) 



In Fig. H] the Chapman-Enskog solution to ( |46T > is plotted to illustrate the relaxation of system 
(OTT i. starting from the following initial pivot and tangent space: 



Y 1 = [0,0, 1.9,0.85f, A = 



1 
1 

-0.276 -1.405 

0.225 0.0282 



(50) 



Finally, the SIM parameterization is chosen as follows: %\ = a, £2 = c 2- We stress that, in 
the computations, the steady state (Y ss , A ss ) does not depend on the initial choice of F 1 and A 1 
(see Fig. [7]). The latter observation is consistent with the idea behind the RRM IT3l . which can 
be elucidated by saying that states on the SIM represent stable steady states of the dynamical 
system (fJJJ). 

The reduced system for the above example (l43i rules the evolution of the slow variables: 



lie 



fi - fi = - c l 

if = h = ~ 2c 2 
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(51) 




Time Time 



Figure 8: (Color online). Comparison between the detailed solution of the system J43I from the state ci = 1.5, cj = 1.5, 
C3 = 2, C4 = (out of SIM) and the reduced solution of the system i5i \ from the state c\ = 1.5, = 1.5. The numerical 
stiff solver odelSs (H] is used with a> = 3, e = 0.025 and r = 3 X 10~ 10 . 



r 


ss 
L l 


Si 

h. 


ss 


ss 
C 4 


i x icr 


13 


-0.6 


-0.85 


0.551271605706663 


0.010023694645548 


1 x 10" 


12 


-0.6 


-0.85 


0.551271602038736 


0.010023695783466 


1 x 10" 


11 


-0.6 


-0.85 


0.551271598194113 


0.010023697018395 
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0.551271605620228 


0.010023694672833 
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-8 











Table 1: The steady state (c", c£ s , c?, cj s ) of 1211 is computed for several choices of the parameter T starting from the 
state ci = -0.6, c% = -0.85, C3 = —1, C4 = 0.5, with ca = 3 and e = 0.025. No convergence is observed for t > 1 X 10 -8 . 



whereas fast variables can be reconstructed by means of the mappings C3 = c^{c\,C2), c\ — 
ca(c\,c 2 ) (Fig. |5l). 

In this respect, in Fig. [8] a solution trajectory of (l43l is compared to the trajectory of (IBTt . 
where the reconstruction of C3 and c 4 is performed using both the exact Chapman-Enskog solu- 
tion (C3 = 0\, C4 = 62) and a linear look up table based on the nodes refined by the linearized 
RRM (see also Figs. [5]and[6]). 

In Table [TJ we report a sensitivity analysis with respect to the parameter r. We notice that an 
estimate of the time scales of a dynamical system can be obtained by a spectral decomposition 
of the Jacobi matrix. For the case in Fig. [5] at equilibrium {c\ — cn — C3 = 0, c\ — 0.25): 

7.45 x 10~ 9 , 
1.49 x 10~ 8 , 



Tl = T2 — 



= -=3.72X10 
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Figure 9: (Color online). Model for the hydrogen-air combustion (3^1 in a closed system under fixed enthalpy (500kJ/kg) 
and pressure (Ibar). Initialization is accomplished using 1371 and {39), where the matrix B is chosen according to the 
spectral quasi equilibrium parameterization Q|. Here, an arbitrary state on a two dimensional QEM is driven on a two 
dimensional SIM by solving the governing equations of the linearized RRM )2U , with n = 9, q = 2, r = 3, r = 1 X 10~' 
using the stiff numerical scheme ode 15s readily available in Matlab® |43| . 



where A; denotes the ;-th eigenvalue. Hence, in the above computations we use t = 3 x 10~ 10 . 
However, the latter parameters was varied within a wide range of values and no significant effect 
was noticed on both the stability of ( f2Tb and the value of its steady state. 

In addition, we test the equations d3"TT ) for computing the mapping d25b describing the local 
fast subspace. To this end, we make use of d40b and d42l) (with D — D — 1, r — 0). 

At steady state of (fJTJ, we observe (at any node of the grid in Fig. [5]) that the columns of the 
matrix A span the subspace defined by the vectors: 

[0,0,1,0], [0,0,0,1], 

in accordance with the assumption that C3 and C4 are the fast variables. 

Finally, the governing equations of the linearized RRM (f2Tb were applied to a more compli- 



cated case of the detailed reaction mechanism for combustion of hydrogen in air 113611 . Here, Eq. 
(|2TT > were tested for computing states of the SIM with dimensions up to q — 5. In Fig. [9] we 
report an example with q -2. 

Moreover, we observed that any steady state of (|3H corresponds to a matrix A whose columns 
are linear combinations of the fast eigenvectors of J ss (i.e., eigenvectors corresponding to the 
n — q — r largest eigenvalues in absolute value). 
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5. Conclusions and outlook 



Based on the rationale behind the film equation <(5j and the Relaxation Redistribution Method 
(RRM), a set of ordinary differential equations (ODEs) is obtained with the aim of mimicking the 
only fast relaxation of a multiscale dynamical system towards a slow invariant manifold (SIM). 
This approach is characterized by a straightforward implementation consisting in solving the 
ODEs by state of the art stiff numerical schemes, and it proves useful for constructing accurate 
approximations of SIM in any dimensions. 

It is worth stressing that, like RRM, convergence of equations (1211) towards a steady state might 
be used for a fully adaptive construction of heterogeneous (i.e. varying dimension in different 
regions of the phase-space) slow invariant manifolds. 

This work sheds light on the connection between the RRM method and the solution of both 
the invariance and film equations as postulated in ||5, 13, 12] (see discussion at the end of section 

ED- 

In addition, the novel algorithm (l3~TT l for approximating the fast subspace is suggested, and a 
possible initialization procedure for both (12H and (|3H is proposed. 

The methods are tested in the case of detailed combustion of hydrogen and air, as well as 
in a benchmark problem of a model with exact Chapman-Enskong solution of the invariance 
equation. 

We stress that, although the presented methodology has been tested in the case of dissipative 
systems with a unique steady state (see Section @), in this paper we show that the governing 
equations (f2Tb and OTI ) of the linearized RRM are based on the general notions of film equa- 
tion (0) and SIM (local) parameterization. Therefore, investigations on the performance of the 
presented method in more general systems with multiple steady states and chaotic behavior (see, 
e -g-> Jl7l I23I0 are planned for future publications. However, in the latter case, new initializa- 
tion procedures are needed since the method discussed in Section [3~4l is suitable for dissipative 
systems equipped with thermodynamic Lyapunov function. 

Finally, it is worth noticingthat t he p roposed approach represents only one possible implemen- 
tation of the RRM method (j5L ll3lll2l1 ). More accurate descriptions of the SIM (to be addressed 
in future publications) can be obtained as well, abandoning the present linear expressions © and 
(fT2T > in favor of higher order approximations. 
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